05 Resumen pronósticos

Author

Eddie Aguilar

library("easypackages")
packages("tidyverse","fpp3", "tsibble", "feasts","fable", "patchwork")
Loading required package: tidyverse
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.0     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.1     ✔ tibble    3.2.0
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Loading required package: fpp3

── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──

✔ tsibble     1.1.3     ✔ fable       0.3.2
✔ tsibbledata 0.4.1     ✔ fabletools  0.3.2
✔ feasts      0.3.0     

── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
✖ lubridate::date()    masks base::date()
✖ dplyr::filter()      masks stats::filter()
✖ tsibble::intersect() masks base::intersect()
✖ tsibble::interval()  masks lubridate::interval()
✖ dplyr::lag()         masks stats::lag()
✖ tsibble::setdiff()   masks base::setdiff()
✖ tsibble::union()     masks base::union()

Loading required package: patchwork

All packages loaded successfully
library("USgas")

Flujo de trabajo limpio de pronóstico

  • Preparación de los datos (limpieza)
  • Gráfica de los datos (visualización)
  • Definición del modelo (especificación)
  • Entrenamiento del modelo (estimación)
  • Revisar el desempeño del modelo (evaluación)
  • Producir pronósticos

Preparación de los datos (limpieza)

Cargar datos en R, o incluso omitir los datos vacíos (NA), filtrado de la serie (mediante paquetería Tsibble).

Gráfica de los datos (visualización)

global_economy %>%
  filter(Country == "Sweden") %>%
  autoplot(GDP) +
    ggtitle("PIB de Suecia") + ylab("$US billions")

Definición del modelo (especificación)

Describir el modelo. Hay distintos tipos, se usará por ejemplo, un modelo lineal de series de tiempo con TLSM:

TSLM(GDP  ~ trend())
<TSLM model definition>

Entrenamiento del modelo (estimación)

Una vez se define el modelo, lo que sigue es entrenar al modelo. Lo que significa pasarle los datos para que, estdísticamente, encuentre los parámetros que realizan el mejor ajuste posible.

fit <- global_economy %>%
  model(Modelo_tendencia = TSLM(GDP ~ trend()))
Warning: 7 errors (1 unique) encountered for Modelo_tendencia
[7] 0 (non-NA) cases

Ahora tenemos un modelo lineal ajustado y el objeto resultante es un mable (model table).

fit
# A mable: 263 x 2
# Key:     Country [263]
   Country             Modelo_tendencia
   <fct>                        <model>
 1 Afghanistan                   <TSLM>
 2 Albania                       <TSLM>
 3 Algeria                       <TSLM>
 4 American Samoa                <TSLM>
 5 Andorra                       <TSLM>
 6 Angola                        <TSLM>
 7 Antigua and Barbuda           <TSLM>
 8 Arab World                    <TSLM>
 9 Argentina                     <TSLM>
10 Armenia                       <TSLM>
# … with 253 more rows

Revisar el desempeño del modelo (evaluación)

Ver como se ajusta el modelo a los datos y compararlo con otros modelos.

Producir pronósticos

Una vez tenemos un modelo ajustado, podemos hacer pronósticos. Usamos forecast(), h = periodo a pronosticar

fcst <- fit %>% fabletools::forecast(h = 36) # h = "3 years"
fcst
# A fable: 9,468 x 5 [1Y]
# Key:     Country, .model [263]
   Country     .model            Year                 GDP        .mean
   <fct>       <chr>            <dbl>              <dist>        <dbl>
 1 Afghanistan Modelo_tendencia  2018 N(1.6e+10, 1.3e+19) 16205101654.
 2 Afghanistan Modelo_tendencia  2019 N(1.7e+10, 1.3e+19) 16511878141.
 3 Afghanistan Modelo_tendencia  2020 N(1.7e+10, 1.3e+19) 16818654627.
 4 Afghanistan Modelo_tendencia  2021 N(1.7e+10, 1.3e+19) 17125431114.
 5 Afghanistan Modelo_tendencia  2022 N(1.7e+10, 1.3e+19) 17432207600.
 6 Afghanistan Modelo_tendencia  2023 N(1.8e+10, 1.3e+19) 17738984086.
 7 Afghanistan Modelo_tendencia  2024 N(1.8e+10, 1.3e+19) 18045760573.
 8 Afghanistan Modelo_tendencia  2025 N(1.8e+10, 1.3e+19) 18352537059.
 9 Afghanistan Modelo_tendencia  2026 N(1.9e+10, 1.3e+19) 18659313546.
10 Afghanistan Modelo_tendencia  2027 N(1.9e+10, 1.3e+19) 18966090032.
# … with 9,458 more rows

Métodos sencillos de pronóstico

Conocidos como benchmark son básicos pero pueden ser muy utiles.

bricks <- aus_production %>% filter_index("1970" ~ "2004")
bricks
# A tsibble: 137 x 7 [1Q]
   Quarter  Beer Tobacco Bricks Cement Electricity   Gas
     <qtr> <dbl>   <dbl>  <dbl>  <dbl>       <dbl> <dbl>
 1 1970 Q1   387    6807    386   1049       12328    12
 2 1970 Q2   357    7612    428   1134       14493    18
 3 1970 Q3   374    7862    434   1229       15664    23
 4 1970 Q4   466    7126    417   1188       13781    20
 5 1971 Q1   410    7255    385   1058       13299    19
 6 1971 Q2   370    8076    433   1209       15230    23
 7 1971 Q3   379    8405    453   1199       16667    28
 8 1971 Q4   487    7974    436   1253       14484    24
 9 1972 Q1   419    6500    399   1070       13838    24
10 1972 Q2   378    7119    461   1282       15919    34
# … with 127 more rows
bricks %>% autoplot(Bricks)

Método del promedio (media)

En este pronósitco, las predicciones de todos los valores futuros son la media de los datos históricos.

bricks %>% model(MEAN(Bricks))
# A mable: 1 x 1
  `MEAN(Bricks)`
         <model>
1         <MEAN>

Método ingenuo (Naive method)

Se toma el último valor como el pronóstico para todos los valores futuros. También se les conoce como pronósitcos de caminata aleatoria, ya que son óptimos en una serie que siga una.

bricks %>% model(NAIVE(Bricks),
                 RW(Bricks)) # hace exactamente lo mismo que NAIVE()
# A mable: 1 x 2
  `NAIVE(Bricks)` `RW(Bricks)`
          <model>      <model>
1         <NAIVE>      <NAIVE>

Método ingenuo estacional (seasonal Naive)

Se le agrega un componente a Naive para lidiar con datos altamente estacionales.

bricks %>% model(SNAIVE(Bricks))
# A mable: 1 x 1
  `SNAIVE(Bricks)`
           <model>
1         <SNAIVE>
vic_elec %>% 
  autoplot(Demand) %>% 
  plotly::ggplotly()
vic_elec %>% 
  model(SNAIVE(Demand)) %>% 
  forecast(h = "1 day") %>% 
  autoplot(vic_elec %>% filter_index("2014-12-30" ~ .))
Warning: There was 1 warning in `filter()`.
ℹ In argument: `time_in(Time, ...)`.
Caused by warning:
! Argument 'roll' is deprecated. Deprecated in version '1.8.4'.

Método de la deriva (drift)

Variación de Naive que permite que aumente o disminuya el pronótico en el tiempo. El aumento del cambio es el cambio promedio en los datos históricos.

bricks %>% model(
  RW(Bricks ~ drift())
)
# A mable: 1 x 1
  `RW(Bricks ~ drift())`
                 <model>
1          <RW w/ drift>

Ejemplos

# Set training data from 1992 to 2006
train <- aus_production %>% filter_index("1992 Q1" ~ "2006 Q4")
# Fit the models
beer_fit <- train %>%
  model(
    Mean             = MEAN(Beer),
    `Naïve`          = NAIVE(Beer),
    `Seasonal naïve` = SNAIVE(Beer),
    Drift            = RW(Beer ~ drift())
  )
# Generate forecasts for 14 quarters
beer_fc <- beer_fit %>% forecast(h= 14)

# Plot forecasts against actual values
beer_fc %>%
  autoplot(filter_index(aus_production, "1992 Q1" ~ .), level = NULL) +
  ggtitle("Forecasts for quarterly beer production") +
  xlab("Year") + ylab("Megalitres") +
  guides(colour=guide_legend(title="Forecast")) +
  geom_vline(xintercept = as.Date("2007-01-01"), color = "firebrick",
             linetype = "dashed") +
  annotate("label", x = c(as.Date("2003-01-01"),as.Date("2009-01-01")),
           y = 550, label = c("Train set", "Test set"),
           color = c("black","blue"))

beer_fc %>%
  filter(.model == "Seasonal naïve") %>% 
  autoplot(filter_index(aus_production, "1992 Q1" ~ .)) +
  ggtitle("Seasonal naive forecast for quarterly beer production") +
  xlab("Quarter") + ylab("Megalitres") +
  geom_vline(xintercept = as.Date("2007-01-01"), color = "firebrick", linetype = "dashed") +
  annotate("label", x = c(as.Date("2003-01-01"), as.Date("2009-01-01")), y = 550, label = c("Train set", "Test set"), color = c("black", "blue"))

Otro ejemplo:

gafa_stock %>% distinct(Symbol)
# A tibble: 4 × 1
  Symbol
  <chr> 
1 AAPL  
2 AMZN  
3 FB    
4 GOOG  
# Re-index based on trading days
google_stock <- gafa_stock %>%
  filter(Symbol == "GOOG") %>%
  mutate(day = row_number()) %>%
  update_tsibble(index = day, regular = TRUE)

# Filter the year of interest
google_2015 <- google_stock %>% filter(year(Date) == 2015)
# Fit the models
google_fit <- google_2015 %>%
  model(
    Mean    = MEAN(Close),
    `Naïve` = NAIVE(Close),
    Drift   = NAIVE(Close ~ drift()),
    SNAIVE  = SNAIVE(Close)
  )
Warning: 1 error encountered for SNAIVE
[1] Non-seasonal model specification provided, use RW() or provide a different lag specification.
google_fit
# A mable: 1 x 5
# Key:     Symbol [1]
  Symbol    Mean   Naïve         Drift       SNAIVE
  <chr>  <model> <model>       <model>      <model>
1 GOOG    <MEAN> <NAIVE> <RW w/ drift> <NULL model>
# Produce forecasts for the 19 trading days in January 2015
google_fc <- google_fit %>% forecast(h = 19)

google_fc
# A fable: 76 x 5 [1]
# Key:     Symbol, .model [4]
   Symbol .model   day        Close .mean
   <chr>  <chr>  <dbl>       <dist> <dbl>
 1 GOOG   Mean     505 N(602, 6766)  602.
 2 GOOG   Mean     506 N(602, 6766)  602.
 3 GOOG   Mean     507 N(602, 6766)  602.
 4 GOOG   Mean     508 N(602, 6766)  602.
 5 GOOG   Mean     509 N(602, 6766)  602.
 6 GOOG   Mean     510 N(602, 6766)  602.
 7 GOOG   Mean     511 N(602, 6766)  602.
 8 GOOG   Mean     512 N(602, 6766)  602.
 9 GOOG   Mean     513 N(602, 6766)  602.
10 GOOG   Mean     514 N(602, 6766)  602.
# … with 66 more rows
# A better way using a tsibble to determine the forecast horizons
google_jan_2016 <- google_stock %>%
  filter(yearmonth(Date) == yearmonth("2016 Jan"))
google_fc <- google_fit %>% forecast(google_jan_2016)
# Plot the forecasts
google_fc %>%
  autoplot(google_2015, level = NULL) +
    autolayer(google_jan_2016, Close, color='black') +
    ggtitle("Google stock (daily ending 31 Dec 2015)") +
    xlab("Day") + ylab("Closing Price (US$)") +
    guides(colour=guide_legend(title="Forecast"))
Warning: Removed 19 rows containing missing values (`()`).

p1 <- google_fc %>%
  filter(.model == "Drift") %>% 
  autoplot(google_2015) +
    autolayer(google_jan_2016, Close, color='black') +
    ggtitle("Google stock (daily ending 31 Dec 2015)") +
    xlab("Day") + ylab("Closing Price (US$)") +
    guides(colour=guide_legend(title="Forecast"))

p2 <- google_fc %>%
  filter(.model == "Naïve") %>% 
  autoplot(google_2015) +
    autolayer(google_jan_2016, Close, color='black') +
    ggtitle("Google stock (daily ending 31 Dec 2015)") +
    xlab("Day") + ylab("Closing Price (US$)") +
    guides(colour=guide_legend(title="Forecast"))

p1 / p2

Valores ajustados (fitted) y residuales

  • Valores reales: Datos históricos de la serie de tiempo.
  • Valores ajustados: Valores (fitted) que tratan de pronosticar los valores reales
  • Residuales: Valores que el modelo no logró capturar. la diferencia entre valores reales y ajustados.

augment() obtiene los valores ajustados y residuales de un modelo:

augment(beer_fit)
# A tsibble: 240 x 6 [1Q]
# Key:       .model [4]
   .model Quarter  Beer .fitted .resid .innov
   <chr>    <qtr> <dbl>   <dbl>  <dbl>  <dbl>
 1 Mean   1992 Q1   443    436.   6.55   6.55
 2 Mean   1992 Q2   410    436. -26.4  -26.4 
 3 Mean   1992 Q3   420    436. -16.4  -16.4 
 4 Mean   1992 Q4   532    436.  95.6   95.6 
 5 Mean   1993 Q1   433    436.  -3.45  -3.45
 6 Mean   1993 Q2   421    436. -15.4  -15.4 
 7 Mean   1993 Q3   410    436. -26.4  -26.4 
 8 Mean   1993 Q4   512    436.  75.6   75.6 
 9 Mean   1994 Q1   449    436.  12.6   12.6 
10 Mean   1994 Q2   381    436. -55.4  -55.4 
# … with 230 more rows

Diagnóstico de residuales

Es importante analizar los residuales de un modelo. Si los residuales presentan patrones, es signo de que el modelo se puede mejorar.

Los residuales de un buen modelo deben de:

  • No estar autocorrelacionadas. Si se detectan correlaciones entre residuos,todavía hay información útil que se debe modelar.
  • La media de los residuos es cero. Si es distinta a cero, el pronóstico está sesgado.

Opcionales: - Los residuos tienen una varianza constante. - Los residuos se distribuyen de manera normal.

Transformaciones Box-Cox pueden ayudar a mejorar los residuales

google_2015 %>% autoplot(Close) +
  xlab("Day") + ylab("Closing Price (US$)") +
  ggtitle("Google Stock in 2015")

aug <- google_2015 %>% 
  model(NAIVE(Close)) %>% 
  augment()

aug
# A tsibble: 252 x 7 [1]
# Key:       Symbol, .model [1]
   Symbol .model         day Close .fitted  .resid  .innov
   <chr>  <chr>        <int> <dbl>   <dbl>   <dbl>   <dbl>
 1 GOOG   NAIVE(Close)   253  522.     NA   NA      NA    
 2 GOOG   NAIVE(Close)   254  511.    522. -10.9   -10.9  
 3 GOOG   NAIVE(Close)   255  499.    511. -11.8   -11.8  
 4 GOOG   NAIVE(Close)   256  498.    499.  -0.855  -0.855
 5 GOOG   NAIVE(Close)   257  500.    498.   1.57    1.57 
 6 GOOG   NAIVE(Close)   258  493.    500.  -6.47   -6.47 
 7 GOOG   NAIVE(Close)   259  490.    493.  -3.60   -3.60 
 8 GOOG   NAIVE(Close)   260  493.    490.   3.61    3.61 
 9 GOOG   NAIVE(Close)   261  498.    493.   4.66    4.66 
10 GOOG   NAIVE(Close)   262  499.    498.   0.915   0.915
# … with 242 more rows
# aug %>% 
#   features(.resid, mean(.,na.rm = TRUE))

aug %>% pull(.resid) %>% mean(na.rm = TRUE) 
[1] 0.9439931
aug %>% autoplot(.resid) + xlab("Día") + ylab("") +
  ggtitle("Residuales del método naïve")
Warning: Removed 1 row containing missing values (`geom_line()`).

Parece que la media es casi cero y la variación es constante, excepto un outlier.

aug %>%
  ggplot(aes(x = .resid)) +
  geom_histogram() +
  ggtitle("Histograma de los residuales")
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

aug %>% ACF(.resid)
# A tsibble: 23 x 4 [1]
# Key:       Symbol, .model [1]
   Symbol .model            lag      acf
   <chr>  <chr>        <cf_lag>    <dbl>
 1 GOOG   NAIVE(Close)        1  0.0976 
 2 GOOG   NAIVE(Close)        2 -0.0726 
 3 GOOG   NAIVE(Close)        3 -0.0748 
 4 GOOG   NAIVE(Close)        4 -0.0433 
 5 GOOG   NAIVE(Close)        5 -0.0398 
 6 GOOG   NAIVE(Close)        6  0.0206 
 7 GOOG   NAIVE(Close)        7 -0.0652 
 8 GOOG   NAIVE(Close)        8 -0.0291 
 9 GOOG   NAIVE(Close)        9 -0.0379 
10 GOOG   NAIVE(Close)       10 -0.00687
# … with 13 more rows
aug %>% ACF(.resid) %>% autoplot() + ggtitle("ACF of residuals")

Se puede ver que los residuos no están autocorrelacionados (ACF), por lo tanto este método (NAIVE) es bueno para esta serie de tiempo, y los pronósticos derivados de este método pueden ser muy buenos.

Todo en una sola gráfica:

google_2015 %>% 
  model(NAIVE(Close)) %>% 
  gg_tsresiduals()
Warning: Removed 1 row containing missing values (`geom_line()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

aus_production %>% 
  filter_index("1992" ~ .) %>%
  model(SNAIVE(Beer)) %>% 
  gg_tsresiduals()
Warning: Removed 4 rows containing missing values (`geom_line()`).
Warning: Removed 4 rows containing missing values (`geom_point()`).
Warning: Removed 4 rows containing non-finite values (`stat_bin()`).

aus_production %>% 
  filter_index("1992" ~ .) %>%
  gg_tsdisplay(Beer)

Veamos con el método de Media:

google_2015 %>% 
  model(MEAN(Close)) %>% 
  augment() %>% 
  pull(.resid) %>% 
  mean(na.rm = TRUE)
[1] 2.710222e-14
google_2015 %>% 
  model(MEAN(Close)) %>% 
  gg_tsresiduals() + 
  ggtitle("Diagnóstico de residuales para el modelo de Media")

Como podemos ver, el método de la media no es óptimo para esta serie de acciones, siendo que la media es 0, al graficar los residuales vemos que muestra un patón (igual al de los datos menos la media), existen claramente autocorrelaciones de una caminata aleatoria y no hay ninguna distribución normal. estos tres factores lo hacen un método nada efectivo para acciones.

Tests de Portmanteau de autocorrelación

Tests para analizar de una forma más formal la presencia o ausencia de autocorrelación en los residuales. Nos ayuda a ver si las primeras \(h\) autocorrelaciones son significativamente distintas de cero o no.

Test de Box-Pierce

Sumatoria de \(r^2_k\) donde \(h\) es el rezago máximo a considerar y \(T\) es la cantidad de observaciones en la muestra. Si \(r^2_k\) es pequeña, entonces el resultado será pequeño. Se recomienda usar \(h = 10\) para datos no estacionales, \(h = 2m\) para datos estaionales (m es el periodo estacional) o como máximo \(h = T/5\).

Hay que tomar en cuenta, que la prueba no es tan buena cuando \(h\) es grande (mayor a \(T/5\)).

Test de Ljung-Box

Generalmente más preciso que Box-Pierce, de igual forma, valores grandes en el resultado indica que no es ruido blanco y los residuales sí están autocorrelacionados.

Hipótesis nula

La hipótesis nula dice que la serie en cuestión no está autocorrelacionada, o sea, es riudo blanco. Para unos residuales, que la hipótesis nula sea cierta es bueno, ya que significa que no están correlacionados los residuales.

¿Cómo saber si es cierta o no?

Siendo \(\alpha\) el nivel de significancia o nivel máximo de error dispuestos a aceptar (usaremos 5%) y p-value un avlor resultante de las pruebas:

Si p-value es menor a \(\alpha\), entonces rechazamos la hiótesis nula, de lo contrario, la aceptamos.

O sea:

  • \(p-value > \alpha\), H0 es cierta, residuales no están autocorrelacionados.
  • \(p-value < \alpha\), H0 es falsa, residuales están autocorrelacionados.

Ejemplos:

# lag=h and fitdf=K
aug %>% features(.resid, box_pierce, lag=10, dof=0)
# A tibble: 1 × 4
  Symbol .model       bp_stat bp_pvalue
  <chr>  <chr>          <dbl>     <dbl>
1 GOOG   NAIVE(Close)    7.74     0.654
aug %>% features(.resid, ljung_box, lag=10, dof=0)
# A tibble: 1 × 4
  Symbol .model       lb_stat lb_pvalue
  <chr>  <chr>          <dbl>     <dbl>
1 GOOG   NAIVE(Close)    7.91     0.637

Como podemos ver para el método Naive, el p-value es mayor al 5%, por lo tanto sus residuos no están autocrrelacionados, son ruido blanco.

google_2015 %>% 
  model(MEAN(Close)) %>% augment() %>% 
  features(.resid, box_pierce, lag = 10, dof = 0)
# A tibble: 1 × 4
  Symbol .model      bp_stat bp_pvalue
  <chr>  <chr>         <dbl>     <dbl>
1 GOOG   MEAN(Close)   2024.         0
google_2015 %>% 
  model(MEAN(Close)) %>% augment() %>% 
  features(.resid, ljung_box, lag = 10, dof = 0)
# A tibble: 1 × 4
  Symbol .model      lb_stat lb_pvalue
  <chr>  <chr>         <dbl>     <dbl>
1 GOOG   MEAN(Close)   2083.         0

Para el caso de Google, no se logra distinguir los residuales del ruido blanco.

Intervalos de predicción

\(c\) es el porcentaje de cobertura de probabilidad. Normalmente se usará 80% y 95%.

Al tener incertidumbre en el ponóstico, entre maoyr sea el intervalo, menos preciso será nuestro pronóstico.

Intervalos de predicción de un paso

Cunado ser realizan pronósitocs de un paso, la desviación estándar del pronóstico es prácticamente la misma que la desviación estándar de los residuos.

Intervalos de predicción de paso múltiple

Entre más aumenta el horizonte de pronóstico, mayor será el intervalo de pronóstico, Esto es, \(\simga_h\) incrementa con \(h\), entonces necesitamos estimaciones de \(\sigma_h\).

Métodos de referencia

Usando fable:

google_2015 %>%
  model(NAIVE(Close)) %>%
  forecast(h = 10) %>%
  hilo()
# A tsibble: 10 x 7 [1]
# Key:       Symbol, .model [1]
   Symbol .model         day        Close .mean                  `80%`
   <chr>  <chr>        <dbl>       <dist> <dbl>                 <hilo>
 1 GOOG   NAIVE(Close)   505  N(759, 125)  759. [744.5400, 773.2200]80
 2 GOOG   NAIVE(Close)   506  N(759, 250)  759. [738.6001, 779.1599]80
 3 GOOG   NAIVE(Close)   507  N(759, 376)  759. [734.0423, 783.7177]80
 4 GOOG   NAIVE(Close)   508  N(759, 501)  759. [730.1999, 787.5601]80
 5 GOOG   NAIVE(Close)   509  N(759, 626)  759. [726.8147, 790.9453]80
 6 GOOG   NAIVE(Close)   510  N(759, 751)  759. [723.7543, 794.0058]80
 7 GOOG   NAIVE(Close)   511  N(759, 876)  759. [720.9399, 796.8202]80
 8 GOOG   NAIVE(Close)   512 N(759, 1002)  759. [718.3203, 799.4397]80
 9 GOOG   NAIVE(Close)   513 N(759, 1127)  759. [715.8599, 801.9001]80
10 GOOG   NAIVE(Close)   514 N(759, 1252)  759. [713.5329, 804.2272]80
# … with 1 more variable: `95%` <hilo>
google_2015 %>%
  model(NAIVE(Close)) %>%
  forecast(h = 10) %>%
  autoplot(google_2015)

Intervalos de predicción con residuales bootstrap

Cuando no es razonable asumir normlidad en los residuos, podemos aplicarles bootstraping, ya que solo asume la no autocorrelación.

Se obteienen muchos escenarios futuros posibles. Para ver algunos, se usa generate:

fit <- google_2015 %>%
  model(NAIVE(Close))

sim <- fit %>%  generate(h = 30, times = 5, bootstrap = TRUE, seed = 123) # only works with dev version

sim
# A tsibble: 150 x 5 [1]
# Key:       Symbol, .model, .rep [5]
   Symbol .model         day .rep   .sim
   <chr>  <chr>        <dbl> <chr> <dbl>
 1 GOOG   NAIVE(Close)   505 1      744.
 2 GOOG   NAIVE(Close)   506 1      747.
 3 GOOG   NAIVE(Close)   507 1      733.
 4 GOOG   NAIVE(Close)   508 1      737.
 5 GOOG   NAIVE(Close)   509 1      739.
 6 GOOG   NAIVE(Close)   510 1      733.
 7 GOOG   NAIVE(Close)   511 1      728.
 8 GOOG   NAIVE(Close)   512 1      730.
 9 GOOG   NAIVE(Close)   513 1      721.
10 GOOG   NAIVE(Close)   514 1      713.
# … with 140 more rows

Generamos 5 escenarios futuros posibles para los siguientes 30 días de trading.

google_2015 %>%
  ggplot(aes(x = day)) +
  geom_line(aes(y = Close), size = 1) +
  geom_line(aes(y = .sim, colour = as.factor(.rep)), data = sim, size = 1) +
  ggtitle("Google closing stock price") +
  guides(col = FALSE)
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
of ggplot2 3.3.4.

Con esto, podemos obtener intervalos de predicción, al calcular los percentiles de los escenarios futuros. El resultado es el ntervalo de predicción bootstrapped. Esto se puede lograr con forecast:

fc <- fit %>% forecast(h = 30, bootstrap = TRUE)
fc
# A fable: 30 x 5 [1]
# Key:     Symbol, .model [1]
   Symbol .model         day        Close .mean
   <chr>  <chr>        <dbl>       <dist> <dbl>
 1 GOOG   NAIVE(Close)   505 sample[5000]  759.
 2 GOOG   NAIVE(Close)   506 sample[5000]  759.
 3 GOOG   NAIVE(Close)   507 sample[5000]  759.
 4 GOOG   NAIVE(Close)   508 sample[5000]  759.
 5 GOOG   NAIVE(Close)   509 sample[5000]  759.
 6 GOOG   NAIVE(Close)   510 sample[5000]  759.
 7 GOOG   NAIVE(Close)   511 sample[5000]  759.
 8 GOOG   NAIVE(Close)   512 sample[5000]  759.
 9 GOOG   NAIVE(Close)   513 sample[5000]  759.
10 GOOG   NAIVE(Close)   514 sample[5000]  759.
# … with 20 more rows
fc %>% autoplot(google_2015) +
  ggtitle("Google closing stock price")

fc <- fit %>% forecast(h = 30)
fc
# A fable: 30 x 5 [1]
# Key:     Symbol, .model [1]
   Symbol .model         day        Close .mean
   <chr>  <chr>        <dbl>       <dist> <dbl>
 1 GOOG   NAIVE(Close)   505  N(759, 125)  759.
 2 GOOG   NAIVE(Close)   506  N(759, 250)  759.
 3 GOOG   NAIVE(Close)   507  N(759, 376)  759.
 4 GOOG   NAIVE(Close)   508  N(759, 501)  759.
 5 GOOG   NAIVE(Close)   509  N(759, 626)  759.
 6 GOOG   NAIVE(Close)   510  N(759, 751)  759.
 7 GOOG   NAIVE(Close)   511  N(759, 876)  759.
 8 GOOG   NAIVE(Close)   512 N(759, 1002)  759.
 9 GOOG   NAIVE(Close)   513 N(759, 1127)  759.
10 GOOG   NAIVE(Close)   514 N(759, 1252)  759.
# … with 20 more rows
fc %>% autoplot(google_2015) +
  ggtitle("Google closing stock price")

Pronósticos con transformaciones

Si se le hace un pronóstico a una serie anteriormente transformada (Box-cox por ejemplo), al terminar el pronóstico, se tiene que hacer una transformación inversa para ver la versión original del pronóstico.

Transformación inversa de Box-Cox:

\(y_t=\)

\({exp(w_t), \text{ si }\lambda = 0}.\)

\((\lambda w_t^\lambda+1)^{1/\lambda}, \text{ en otro caso}.\)

La paquetería fable realiza la transofrmación ivnersa en automático, cuando se especifica en la definición del modelo.

Ajustes por sesgo

Un problema al realizar transformaciones matemáticas, como Box–Cox es que las estimaciones puntuales re transformadas ya no presentan la media de la distribución de predicción, sino que la mediana. Esto puede llegar a ser problema dependiendo del contexto, pero normalmente no causa mucha diferencia.

Ajustar por sesgo es la transformación inversa de la media, entre más grande sea la varianza, mayor será la diferencia entre la media y la mediana.

eggs <- as_tsibble(fma::eggs)
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 
eggs %>% 
  model(RW(log(value) ~ drift())) %>% 
  forecast(h=50) %>% 
  autoplot(eggs, level = 80,  
           point_forecast = lst(mean, median))
Warning in ggplot2::geom_point(mapping = mapping, data =
dplyr::semi_join(object, : Ignoring unknown aesthetics: linetype

Por lo tanto, para tener un pronóstico sesgado, se usa point_forecast = lst(median)

Pronósticos con descomposición

La descomposición de series de tiempo puede ser útil para producir pronósticos.

El pronóstico se realiza en dos pasos: - Un pronóstico para el componente estacional. - Un pronósitco separado para la serie desestacionalizada.

En realidad, SNAIVE es el pronóstido de solo el componente estacional.

(us_retail_employment <- fpp3::us_employment %>%
  filter(year(Month) >= 1990, Title == "Retail Trade"))
# A tsibble: 357 x 4 [1M]
# Key:       Series_ID [1]
      Month Series_ID     Title        Employed
      <mth> <chr>         <chr>           <dbl>
 1 1990 Jan CEU4200000001 Retail Trade   13256.
 2 1990 Feb CEU4200000001 Retail Trade   12966.
 3 1990 Mar CEU4200000001 Retail Trade   12938.
 4 1990 Apr CEU4200000001 Retail Trade   13012.
 5 1990 May CEU4200000001 Retail Trade   13108.
 6 1990 Jun CEU4200000001 Retail Trade   13183.
 7 1990 Jul CEU4200000001 Retail Trade   13170.
 8 1990 Aug CEU4200000001 Retail Trade   13160.
 9 1990 Sep CEU4200000001 Retail Trade   13113.
10 1990 Oct CEU4200000001 Retail Trade   13185.
# … with 347 more rows
us_retail_employment %>% 
autoplot(Employed)

dcmp <- us_retail_employment %>%
  model(STL(Employed ~ trend(window = 7), robust=TRUE)) %>%
  components() %>%
  select(-.model)

dcmp
# A tsibble: 357 x 7 [1M]
# Key:       Series_ID [1]
   Series_ID        Month Employed  trend season_year remainder season_adjust
   <chr>            <mth>    <dbl>  <dbl>       <dbl>     <dbl>         <dbl>
 1 CEU4200000001 1990 Jan   13256. 13326.       -72.6     2.85         13328.
 2 CEU4200000001 1990 Feb   12966. 13297.      -273.    -58.1          13239.
 3 CEU4200000001 1990 Mar   12938. 13270.      -295.    -36.3          13233.
 4 CEU4200000001 1990 Apr   13012. 13230.      -216.     -1.59         13228.
 5 CEU4200000001 1990 May   13108. 13223.      -119.      4.72         13228.
 6 CEU4200000001 1990 Jun   13183. 13212.       -30.6     0.978        13213.
 7 CEU4200000001 1990 Jul   13170. 13198.       -32.8     5.23         13203.
 8 CEU4200000001 1990 Aug   13160. 13178.       -20.1     1.84         13180.
 9 CEU4200000001 1990 Sep   13113. 13149.       -45.5     9.76         13159.
10 CEU4200000001 1990 Oct   13185. 13113.        64.2     8.28         13121.
# … with 347 more rows
dcmp %>%
  model(NAIVE(season_adjust)) %>%
  forecast() %>%
  autoplot(dcmp) + ylab("New orders index") +
  ggtitle("Pronóstico naïve de los datos desestacionalizados")

Le podemos agregar nuevamente la estacionalidad con la función decomposition_model():

us_retail_employment %>%
  model(stlf = decomposition_model(
             STL(Employed ~ trend(window = 7), robust = TRUE),
             NAIVE(season_adjust)
  )) %>%
  forecast() %>%
  autoplot(us_retail_employment)

Evaluación del desempeño de los pronósticos

Conjuntos de entrenamiento y prueba

El tamaño de la prueba es, generalmente, del 20% del total de datos disponibles.

  • Un modelo que se ajusta muy bien a los datos de entrenamiento no necesariamente produce los mejore spronósticos.
  • Al aumentar los parámetros, podemos llegar a un ajuste perfecto del modelo a los datos.
  • Puede darse un efecto de sobre ajuste (over-fitting) y esto es tan malo como tener un muy mal ajuste.

Funciones para sementar las series de tiempo

top_n nos permite obtener las n observaciones más extremas.

gafa_stock %>%
  group_by(Symbol) %>%
  top_n(1, Close)
# A tsibble: 4 x 8 [!]
# Key:       Symbol [4]
# Groups:    Symbol [4]
  Symbol Date        Open  High   Low Close Adj_Close   Volume
  <chr>  <date>     <dbl> <dbl> <dbl> <dbl>     <dbl>    <dbl>
1 AAPL   2018-10-03  230.  233.  230.  232.      230. 28654800
2 AMZN   2018-09-04 2026. 2050. 2013  2040.     2040.  5721100
3 FB     2018-07-25  216.  219.  214.  218.      218. 58954200
4 GOOG   2018-07-26 1251  1270. 1249. 1268.     1268.  2405600

Errores de pronóstico

Diferencia entre el valor real ocurrido y el dato pronosticado.

Los errores de pronóstico son distintos de los residuales en:

  • Los residuales se calculan con los datos de entrenamiento, mientras que los errores de pronóstico se calculan con los de prueba.
  • Los residuos se calculan mediante pronósticos de un paso, los errores pueden ser multi-step.

Errores dependientes de la escala de los datos

Dependen de la escala de la serie de tiempo, por lo que no se puede comparar con series de tiempo con otras unidades.

Los dos más utilizados son:

  • MAE, Mean absolute error, el promedio del error.
  • RMSE, Root mean squared error, la raíz cuadrada del promedio del error.

MAE es muy utilizado por su facilidad de cómputo e interpretación.

Errores porcentuales

Son un porcentaje y son muy utilizados para comparar el desempeño de pronósticos de distintos conjuntos de datos.

-MAPE, Mean absolute percentage error, el promedio del error.

El problema es que si un dato es 0, este error será infinito, por lo tanto se inventó:

-sMAPE, symmetric MAPE, resuelve el problema del MAPE.

Sin embargo, no es recomendable en la práctica.

Errores escalados

ALternativa para porcentuales para comparar con otras series.

  • MASE, Mean absolute scalated error, promedio de los errores escalados.
recent_production <- aus_production %>% filter(year(Quarter) >= 1992)
beer_train <- recent_production %>% filter(year(Quarter) <= 2007)

beer_fit <- beer_train %>%
  model(
    Mean = MEAN(Beer),
    `Naïve` = NAIVE(Beer),
    `Seasonal naïve` = SNAIVE(Beer),
    Drift = RW(Beer ~ drift())
  )

beer_fc <- beer_fit %>%
  forecast(h = 10)

beer_fc %>%
  autoplot(recent_production, level = NULL) +
  xlab("Year") + ylab("Megalitres") +
  ggtitle("Forecasts for quarterly beer production") +
  guides(colour=guide_legend(title="Forecast"))

Usamos accuracy para ver los errores en el ajuste de modelos (mabble) a los datos de entrenamiento:

accuracy(beer_fit)
# A tibble: 4 × 10
  .model         .type           ME  RMSE   MAE    MPE  MAPE  MASE RMSSE   ACF1
  <chr>          <chr>        <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>
1 Mean           Training  0         43.6  35.2 -0.937  7.89  2.46  2.60 -0.109
2 Naïve          Training  4.76e- 1  65.3  54.7 -0.916 12.2   3.83  3.89 -0.241
3 Seasonal naïve Training -2.13e+ 0  16.8  14.3 -0.554  3.31  1     1    -0.288
4 Drift          Training -5.41e-15  65.3  54.8 -1.03  12.2   3.83  3.89 -0.241

Ahora la usamos para ver los errores de pronóstico (fabble) comparado a los datos de prueba:

beer_fc |>  accuracy(recent_production)
# A tibble: 4 × 10
  .model         .type    ME  RMSE   MAE    MPE  MAPE  MASE RMSSE    ACF1
  <chr>          <chr> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>   <dbl>
1 Drift          Test  -54.0  64.9  58.9 -13.6  14.6  4.12  3.87  -0.0741
2 Mean           Test  -13.8  38.4  34.8  -3.97  8.28 2.44  2.29  -0.0691
3 Naïve          Test  -51.4  62.7  57.4 -13.0  14.2  4.01  3.74  -0.0691
4 Seasonal naïve Test    5.2  14.3  13.4   1.15  3.17 0.937 0.853  0.132 

Los errores de ajuste (datos de entrenamiento) sirven para ver que modelos utilizar para pronosticar.

Los errores de pronóstico (datos de prueba) nos da mayor claridad de cuáles modelos parecen producir los mejores pronósticos.

Después de esto se recalculan los modelos utilizando toda la información (sin separar en conjuntos de entrenamiento y prueba) y se producen los pronósticos reales a presentar.